from scipy import *
from scipy.integrate import quad
from numpy import *
#import numpy
#from matplotlib import pyplot
k = 2*pi
# contour parametrization and its derivatives [alpha,betta]->R
alpha = 0.0
betta = 1.0
def ro1(t):    return t
def z1(t):    return 0.5*(1-t**2)
def dro1(t):    return t-t+1
def dz1(t):    return -t
def ddro1(t):    return t-t
def ddz1(t):    return t-t
def dddro1(t):    return t-t
def dddz1(t):    return t-t

#discretization order
length = quad(lambda t: numpy.sqrt(dro1(t)**2+dz1(t)**2),alpha,betta)
n = int(3*length[0]*k/numpy.pi) + 5
print "n =",n

# the  roots of Chebyshev polynomials(Chebyshev points)
ind1 = arange(n-1)
ind2 = arange(n)
t0 = numpy.cos(numpy.pi*(ind1+1)/n)
t = numpy.cos(numpy.pi*(2*ind2+1)/(2*n))
ind = arange(2*n-1)
t0t = hstack((t0,t))

# values of contour parametrization 
# and its derivatives [-1,1]->R in the Chebyshev points
mult = (betta-alpha)/2
tau = alpha*(1-t0t)/2 + betta*(1+t0t)/2
ro = ro1(tau)
z = z1(tau)
dz = dz1(tau)*mult
dro = dro1(tau)*mult
ddz = ddz1(tau)*mult**2
ddro = ddro1(tau)*(mult**2)
dddro = dddro1(tau)*(mult**3)
dddz = dddz1(tau)*(mult**3)

#Lyame koefficient
h_tau = numpy.sqrt(dro**2+dz**2)

#distance between to point on the surface
def L(m,p,psi): return numpy.sqrt(ro[p]**2+ro[m]**2-2*ro[p]*ro[m]*cos(psi)+(z[p]-z[m])**2)

#LIMITS
  #Part12 = numpy.exp(-1j*k*LL)/(LL**2)*(-1/LL - 1j*k)*numpy.cos(M*psi)
  #Part22=numpy.exp(-1j*k*LL)/(LL**3)*(-k**2 + 3j*k/LL + 3/(LL**2))*numpy.cos(M*psi)
def f_Smooth_Re(M,p,psi):
    LL = numpy.sqrt(2)*ro[p]*numpy.sqrt(1-cos(psi))
    Part11 = -(dro[p]**2)*numpy.cos(psi)-dz[p]**2
    Part12 = -((numpy.cos(k*LL))/(LL**3) + k*numpy.sin(k*LL)/(LL**2))*numpy.cos(M*psi)
    Part13 = (1/LL**3)*(1 - (M**2 * numpy.sin(psi)**2)/2)
    Part14 = (k**2/2 + M**2/(2*ro[p]**2)*(numpy.cos(psi) - 1))/LL
    Part21 = (ro[p]*dro[p]*(1 - numpy.cos(psi)))**2
    Part22 = (numpy.cos(k*LL)/(LL**3)*(-k**2 + 3/(LL**2)) + 3*k*sin(k*LL)/(LL**4))*numpy.cos(M*psi)
    Part23 = -3*(1 - (M**2 * numpy.sin(psi)**2)/2)/(LL**5)
    Part24 = -(1/LL**3)*(k**2/2 + M**2/(2*ro[p]**2)*(cos(psi)-1))
    Res = Part11*(Part12 + Part13 + Part14) + Part21*(Part22+Part23+Part24)
    return Res
def f_Smooth_Im(M,p,psi):
    LL = numpy.sqrt(2)*ro[p]*numpy.sqrt(1-cos(psi))
    Part11 = -(dro[p]**2)*numpy.cos(psi)-dz[p]**2
    Part12 = ((numpy.sin(k*LL))/(LL**3) - k*numpy.cos(k*LL)/(LL**2))*numpy.cos(M*psi)
    Part21 = (ro[p]*dro[p]*(1 - numpy.cos(psi)))**2
    Part22 = (-numpy.sin(k*LL)/(LL**3)*(-k**2 + 3/(LL**2)) + 3*k*cos(k*LL)/(LL**4))*numpy.cos(M*psi)
    Res = Part11*Part12 + Part21*Part22
    return Res

print(f_Smooth_Im(2,5,1))
def Smooth(M,p):   
    #Int = quad(lambda psi: f_Smooth(2,5,psi),0,numpy.pi)
    Int,err = quad(f_Smooth_Re,0,numpy.pi, args=(2,5,))
    return 2*Int
print Smooth(2,5)
